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PILOT  EFFORT  TO  DEFINE  IMPROVED  NUMERICAL  ANALYSIS  METHODS  FOR  THE 
HANDLING  OF  MODERATE  APERTURE  ACOUSTIC  ARRAY  DATA 


Background 

The  Code  5120  research  project  "High  Resolution  of  Complex  Acoustic 
Fields  Using  Moderate  Apertures"  will,  if  successful,  enhance  the  opera¬ 
tional  Navy's  ability  to  detect  and  track  submarines  in  shallow  water  or 
Arctic  environments •  Initially  the  project  is  considering  the  following 
physical  circumstances:  high  signal  to  noise  ratio;  shallow  water  (300 
meters  maximum);  a  short  array  of  hydrophones  (less  than  300  meters);  and 
a  small  number  of  hydrophones  (up  to  64). 

As  part  of  this  initial  effort  numerical  analysis  and  processing 
techniques  were  to  be  developed.  These  methods  may  then  be  incorporated 
into  experiments  as  part  of  the  signal  processing  portion  of  the  experi¬ 
mental  apparatus.  Thus,  successful  research  should  provide  numerical- 
experimental  methods  to  better  represent  an  incident  sound  pressure  field 
everywhere  along  an  array  of  hydrophones,  the  major  source  of  that  field 
initially  being  controlled  by  the  experimenter. 

(Underlying  the  numerical  analysis  efforts  is  the  assumption  that  we 
will  ultimately  be  involved  in  a  "bootstrapping"  operation.  That  is,  a 
mathematical  model  derived  from  the  experimental  data  will  be  useful  for 
improving  the  collection  and  reduction  of  experimental  data;  in  turn  these 
new  data  serve  as  the  basis  for  further  analysis  and  refinement  of  the 

mathematical  model,  and  so  on.) 

Manuscript  approved  September  20,  1962. 
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The  experiments  will  probably  Involve  placing  sound  sources  with 
specifically  designed  sonic  characteristics  in  the  water  at  preselected 
points.  The  mathematics  utilized  in  the  signal  processing  equipment  depends 
upon  the  results  of  the  research  in  numerical  processing  techniques. 

In  submarine  detection  and  tracking  successful  techniques  would  also 
have  to  deal  with  a  low  signal-to-noise  ratio.  It  is  hoped  the  methods 
developed  in  this  project  will  shed  light  on  techniques  for  use  in  such 
noisy  environments. 

To  initiate  the  development  of  the  numerical  and  processing  techniques) 
a  short  term  (approximately  2  month)  investigation  was  undertaken  by  Charles 
Osgood  in  collaboration  with  Richard  Heitmeyer  and  William  Moseley. 

This  investigation  made  use  of  the  following  assumptions: 

o  The  sound  pressure  field  f(z)  can  he  well  represented  by  a  linear 
combination  of  a  few  modal  solutions  to  the  wave  equation; 

o  The  modal  solutions  are  known  well  numerically  (even  if  not 
analyt ically); 

o  The  eigenvalues  kn  are  also  known  well  numerically; 

o  Each  pressure  field  is  well  approximated  by  a  linear  combination  of, 
say,  the  100  lowest  modals; 

o  Only  approximately  20  modals  are  ever  neelod  to  represent  any  particu¬ 
lar  pressure  field;  and, 

o  The  slgnal-to-noise  ratio  is  high,  but  some  noise  is  always  present. 
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Introduction 


This  paper  presents  for  consideration  a  number  of  approaches  for 
attacking  the  problem  above  which  were  developed  over  a  ten  week  period.  It 
is  hoped  that  further  research  in  this  area  will  be  aided  by  the  Ideas 
presented. 

A  few  words  about  the  organization  of  this  paper:  Section  I  addresses 
the  problem  of  determining  those  modals  un(z)  actually  present  in  the 
representation  of  the  incident  sound  pressure  field  f(z).  Section  II  treats 
the  problem  of  determining  the  coefficients  of  these  modal  functions  u^z). 
There  are  five  appendices  expanding  on  different  points.  The  levels  of 
the  different  appendices  vary  from  being  more  tutorial  than  the  main  text 
to  being  more  technical  than  the  main  text. 


Section  1 


Let  f(z)  denote  the  incident  preseure  field  on  the  water  column  where 
measurements  are  being  made.  Let  H  denote  the  depth  of  the  water.  Let 
the  zj,  1  £  j  £  N,  denote  the  assumed  depths  of  the  N  hydrophones.  As  it 
happens,  we  shall  always  assume  that  the  hydrophones  are  equally  spaced,  so 
let  h  -  Zj+l  -  zj*  It  is  assumed  that  the  emitter  has  a  constant  frequency 
u.  Let  c(z)  denote  the  speed  of  sound  at  depth  z  in  a  neighborhood  of  the 
water  column. 

From  [2]  we  know  that  in  shallow  water  we  may  assume  that 

f(z)  -  £n  An  un(z) 

where  each  modal  function  satisfies  a  differential  equation  of  the  form 


0) 


C2(z) 


u  (z)  *  0 
n 


for  some  constant  k„. 

Determining  the  kjj's  of  the  modals  Involved  in  representing  a  particular 

a2  2 

a  to 

f(z)  is  the  problem  now  considered.  LetUi* — y  “  —5 - •  It  is  shown  in 

dz^  c  (z) 

Appendix  I  that  if  d  modals  are  needed  to  represent  f(z)  then  for  almost  all 
choices  of  equally  spaced  zi,...,z<]  the  polynomial  equation  in  y 

1  f (z^)  ...  f(zd) 

yij>f(zi)  •  •  •  f  ( *d) 


-  0  (1) 


]  yd  ij^fCzp...  ^df(*d)  j 

will  not  vanish  Identically  and  will  have  as  roots  exactly  those  (~k*> 
corresponding  to  the  d  modals  u^z)  needed  to  represent  f(z).  (If  fewer 
than  d  modals  are  needed  then  the  polynomial  in  (1)  vanishes  identically.) 


Obviously  in  trying  to  use  (1)  there  is  (probably  well  founded)  concern 
about  approximating  high  order  derivatives  numerically.  One  should  calculate 
the  polynomial  in  (1)  with  t  replacing  d  for  t  -  1,2, • ••,  until  a  zero 
polynomial  is  obtained.  Then  backing  up  to  a  d+1  by  d+1  matrix,  the  roots 

fl 

of  (1)  should  be  determined.  Error  in  approximating  the  ,])  f(zj)  could  even 

make  it  impossible  to  tell  which  polynomial  equation  to  solve.  (On  the 

% 

other  hand,  if  the  values  of  the  ill  f(zj)  have  been  accurately  determined 
"solved”  is  too  strong  a  word  to  describe  the  necessary  calculations.  The 
roots  must  he  included  among  a  set  of  no  more  than  100  numbers;  therefore, 
the  polynomial  need  only  be  evaluated  100  times.) 

A  variant  of  the  above  idea  is  discussed  in  Appendix  II.  It  seems  more 
likely  that  differences  of  the  form  f(z+h)-f(z)  can  be  more  accurately 
obtained  than  can  derivatives  of  f(z).  The  difficulty  with  this  observation 
is  that  we  seem  to  have  only  a  differential  equation  available  for  use.  In 
a  very  important  case,  however,  we  do  have  difference  equations  which  are 
satisfied  by  the  modal  functions. 

In  the  case  of  constant  sound  speed  (in  the  water  column)  we  obtain 


solutions  of  the  form  un(z)  -  sin(6nz).  [More  correctly  un(z)  is  propor¬ 
tional  to  8ln(8|)Z)  but  we  shall  Ignore  that  point  -in  this  paper.]  Second 
it  is  "folklore"  in  the  field  that  in  cases  of  actual  Interest  the  higher 
subscripted  models  look  as  if  they  are  sin's.  In  Appendix  III  it  is  shown 
that  if  one  replaces  throughout  equation  (1)  by  the  operator  where 


^hS(z)  "  g(z+h)-2g(z)+g(z-h)  then  the  polynomial  equation  replacing  (1)  has 

1  2 

roots  equal  to  the  numbers  -  (2  sin  (8njh))  .  [Since  we  do  not  divide  by 

,  -2  ■  d2 

hz  the  noise  is  less  likely  to  be  exaggerated.  Note  that  h  — *r 

dz 

as  h-*-0  and  that  -h“^(2  sinS^h)*-*-  0fl  as  h  -*-0.] 
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What  abouc  the  lower  subscripted  modals,  in  the  situation  where  the 
higher  subscripted  modals  look  sinusoidal?  Let  un(z)  be  one  of  the  higher 


subscripted  modals  which  is  supposed  to  look  like  sin  @nz  for  some  0n.  Since 

d2  , 

7  sin  0_z  ■  -  0„  sin  0_z  and 

dz  n  n  n 


d 

dz 


"“<z)  (t77T  ' k" )  “»<z> 


we  conclude  that 


0) 


c2(z) 


2  2 

-  k  a  -  0  . 
n  n 


This  approximation  is  probably  less  good  than  the  original  approximation 

un(z)ssin  0nz.  We,  therefore,  draw  only  the  weak  (and  heuristic)  conclusion 

from  this  that  c(z)  varies  slowly.  It  should  then  follow  that  the  lower 

subscripted  modals  look,  locally,  like  sin(<j>n  +0nz),  for  constants  0n  and 

,j)n.  [The  general  solution  of  y"  -  -  02y  is  of  the  form  An  8in(<f>n  +®„z)]* 

Therefore,  if  the  hydrophones  are  on  a  small  enough  interval,  we  would  expect 

1  o 

the  above  polynomial  equation  to  have  roots  -  (2  sin  ©njT  h))  corresponding 
to  lower  subscripted  modals  which  are  (locally)  proportional  to  sin  (q>n  +0nz). 

That  is,  we  would  expect  to  obtain  these  numbers  as  roots  except  for 
noise  in  the  measured  values  of  f(z).  Ideally  we  would  like  to  smooth  the 
da.:  a  rather  than  differencing  it.  In  a  sense  it  is  possible  to  accomplish 
our  desire  by  using  some  mathematical  trickery  applied  to  the  above 
(difference)  method  of  calculating  the  roots  -(2  sin  (9n|h))2. 

In  Appendix  III  it  is  shown  how  we  may  construct  Fjc(z) ,  a  linear  combi¬ 
nation  of  the  functions  f(z),  f(z-h),  f(z-2h),...,  with  positive  coefficients, 
which  when  differenced  still  stays  a  linear  combination  of  the  f(z),  f(z-h), 

f (z-2h) .  with  positive  coefficients.  Further,  one  can  use  the  set  of 

Fk(z)'s  (there  is  one  for  k  ■  0,1,...)  to  construct  a  polynomial  equation 
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from  Whose  roots  the  0n  can  be  determined.  (On  the  other  hand,  this  method 
does  require  more  date  points  then  do  the  methods  discussed  previously; 
additionally,  the  degree  of  the  polynomial  equation  doubles  from  d  to  2d). 

Trying  to  verify  the  locations  of  roots  to  these  polynomial  equations 
may  prove  difficult.  A  probably  much  better  method  of  root  verification 
applies  to  the  polynomial  equation  obtained  in  Appendix  III  ((2)  below)  as 
well  as  to  a  modified  version  of  equation  (1)  ((3)  below);  this  method  Is 
discussed  next. 

The  equation  obtained  in  Appendix  III  Is  of  the  form 
1  fo<*a>  Fl(^>  F2d-l(*£> 

y  f1<*4)  F2<*4>  F2d(*A> 

-  0  (2) 

see  i 

see  « 

y2d  r2i^xO  F2d+l(*i ^  F4d-l(*fc) 

One  could  have  obtained  an  analogue  of  equation  (1)  of  the  form 
1  f<sd)  *f(xd)  ...  *  d"1f(*d) 

y  4/f(xd)  4>2f(*d)  •••  *df(*d) 

.  .  .  •  ■  0  (3) 

me  a  e 

e  a  a  e 

yd  /f(zd)  i ...  4>2d"1  f(td) 

which  has  the  same  roots  as  equation  (1).  [If  4*h  substitutes  for  ty,  (3)  has 

1  9 

the  roota  -  (2  sin  (d^h))*  If  the  modals  are  sinusoidal.  Note  (3)  needs 
more  data  points  than  does  (1).] 

Form  the  polynomials 
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-Vflfr 


■  -  -qwrr^art  jm*-’?-’ *i  r*t  "tt. 


respectively.  Consider  the  problem  of  approximating  these  respective 
polynomials  by  rational  functions.  There  always  exist  polynomials  uj.,  V]_,  * 

and  vj  of  degrees  not  exceeding  d-1,  d,  2d-l,  and  2d,  respectively,  such 
that: 

~1  ^d-l  .  . 

u. v •  =  I  (tyJ f (Zj) )x^ 

1  1  j-0  d 

plus  terms  of  higher  degree  in  x 
and 

-1  4d_1  i 

u,v,  =  E  F  (z.)xJ 

22  j-0  j 

plus  terms  of  higher  degree  in  x.  Almost  always  (in  a  sense  which  could  be 
made  mathematically  precise)  u^,  vj,  U2»  and  V2  are  of  degrees  exactly  d-1, 
d,  2d-l,  and  2d.  Further  these  polynomials  are  unique  (if  we  demand  their 
leading  coefficient  be  one)  and  the  zeros  of  vi  and  V2  are,  respectively,  the 
reciprocals  of  the  roots  of  (3)  and  the  reciprocals  of  the  rootB  of  (2). 

As  is  explained  in  Appendix  IV  there  is  an  efficient  algorithm  to 
calculate  values  of  the  rational  functions  Uj  Vj”*  and  U2  V2~*.  If  we  are 
at  a  pole  (infinity)  of  the  rational  function  we  have  determined  a  root  of 
the  associated  equation. 

The  residues  of  these  rational  functions  at  their  poles  can  be  related 

to  the  coefficients  of  the  models  in  their  representation  of  f(z).  For 

-1  2  -2  lln  2 
example,  the  residue  of  u^  Vj  at  z  ■  He*  is  ~\^-n  •  Therefore  z  -kn 

—  2  —1  —2 

(z  +  kn  )  “i  vi  ■  "An  kjj  .  It  might  or  might  not  be  reasonable  to 
compute  the  A,,  this  way.  In  any  event  this  brings  us  to  consideration  of 
ways  of  calculating  the  An. 
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Section  II 


Finding  the  Aq's  can  be  approached  in  several  ways  depending  upon  what 
assumptions  are  made.  In  the  general  case  we  have  that 


U  (k;  *  k"))va<‘> 


(4) 


[For  a  more  gradual  introduction  to  the  use  of  products  of  differential 
operators  see  Appendix  V.] 

Because  of  the  difficulty  of  computing  derivatives  numerically,  this  may 
not  be  a  feasible  way  of  finding  the  A„.  As  in  Section  I,  it  seems  reasonable 
to  utilize  the  operator  ip),.  Then  we  can  write 

(j!»  (*h+  (2  sin  ej  f  )2))  f<zd-i>  • 

(j*>  ((2sl"  6j!  )2  -  (2sl"  9nl  )))  Vn(za-1>- 


This  expression  should  be  valid  whenever  each  i>n(z)  a  sin(<£n  +0nz).  (Some¬ 
thing  like  this  idea  should  work  in  the  situation  of  Appendix  III  where  f(z) 
is  replaced  by  a  linear  combination  of  f(z),  f(z-h),...  having  positive 
coefficients.  The  exact  algebra  would  be  different  however.) 

There  is  another  approach,  however,  which  is  what  will  be  discussed  next. 
[In  Appendix  V  more  detail  is  provided.]  Set  x  -  ir  (z-zj)(zjj-zj)-1  and  let  k 
denote  a  non-negative  integer.  Suppose  one  multiplies  both  sides  of  (4)  by 
(sin  it x)  u  un(z)  and  integrates  the  resultant  expression  over  [z^,zN]. 
Because  of  the  high  order  of  vanishing  of  the  integrands  at  z  ■  zj  and  z  -  zfj 
it  is  possible  to  use  integration  by  parts  2(d-l)  times  on  the  left  hand 
Integral  and,  finally,  obtain  an  Integrand  there  which  is  free  of  derivatives 
of  f(z).  This  integral  looks  like 
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At  least  the  2d-2  powe”  of  sin  irx  is  mandated  by  the  requirement  that 
we  perform  2d-2  Integrations  by  parts  without  having  derivatives  of  f(z) 
appear.  The  reason  for  possibly  wishing  that  k  >  0  is  that  then  the 
Integrand  in  (5)  will  agree  (in  value  and  derivatives)  up  to  the  order  k 
at  z  *  zj  and  at  z  -  zN.  Since  (5)  must  be  evaluated  by  numerical  inte¬ 
gration  we  are  interested  in  the  rate  of  convergence  of  some  types  of 
numerical  integration  schemes.  Having  the  integrand  agree  to  a  high  order 
at  the  two  points  z  ■  and  z  ■  zN  can  be  expected  to  add  considerable 
accuracy  in  the  numerical  integration.  (See  Appendix  V  for  details.) 

Appendix  V  discusses  approaches  to  evaluating 

[>"-  (t1 "  7^  +kj)]  '”ln  *>2d~2+k 

quickly  using  the  Fast  Fourier  Transform. 


Appendix  I 


A  possible  way  to  determine  the  kn's  corresponding  to  the  nodal  function 
un(z)  involved  in  f(z)  is  the  following: 


Write 


f<*>  ■  Zn  u„(z) 


where  the  ^(z)  are  the  d  modal  functions  appearing  in  the  decomposition  of 
f(z).  Suppose  that  one  can  determine  good  numerical  approximations  to 

f(z)'(kz --$&))  f(z)'(i^-^ur)  f(z) . (^-^%y)df(z> 

for  the  d  different  points  zlt  *2  •  *1  +  h,...,  z<j  ■  +  (d-l)h.  At  each 

point  zj  the  values  of  the  different  /  j 2  \  i  are  of  the  form 

y  2  1  Idz^  c^(z)  J  ^Zj^ 

En  An  (^n)1  un<2j)‘  '  ' 

For  simplicity  resubscript  the  d  modals  so  they  are  denoted  as  ui(z),..., 

o 

u^Cz).  Since  the  -k£  are  distinct  numbers  the  determinant  of  the  matrix 


(He?)*’1  (-kf)** 


(-ki)^1 


does  not  vanish.  (This  is  by  the  nonvanishing  of  the  Vandermonde  determinant.) 


A1  U1  (*l)  A1  “1  <*2>  •••  A1  U1  (*d> 


a2  u2  (*1>  a2  u2  <*2) 


a2  “2  (^d) 


■  <An  un  (*J»- 


Ad  u<j  (zj)  Aj  ty  (Z2) 


Aj  uy  (zj) 


If  B  is  nonsingular  Chen  so  is  the  matrix 


AB 


where  1  <  i  +  1  <  d  and  1  £  j  <  d. 

Since  by  hypthesis  +  0,  the  nonsingularity  of  B  is  equivalent 

to  the  nonsingularity  of  the  matrix  (u^Zj)).  Assume  the  ^  are  analytic 
functions.  Then  lu^zi  +  (j-1)  h)  I  ■  M(z^,h)  is  an  analytic  function  of  h. 
For  each  fixed  zj,  if  M(z^,h)  vanishes  for  an  infinite  number  of  h,  0 
£  H,  then  M(z^,h)  =  0,  for  this  value  of  z^,  as  a  function  of  h.  Define 

Ah  fey  Ah  8<z>  "  (8(*+h)  -  gC*))**”1. 

There  exist  constants  ej£,  depending  upon  h  but  independent  of  n  such 

that: 


/  u1(z1)Ah  ux  (z1)  ...  AJ’1  ux  (zt)  \ 

/  ®11*  * *eld 

d-1 

u2(zl>Ah  u2  (zl>  *•*  A  u2  (*1> 

e  • 

•  •  h  • 

e  e  e 

e  e  e 

1  ud^2^  A h  ud  (*i)  •••  A^  1  ud  (zi^  j 

•  • 

e  e 

\  edl*  *  *edd  t 

/  UX  (Zj)  ...  U1(zd) 


\  ud  (zx)  ...  U,j(zd)  I 
hence,  M(zi,h)  -  N(h)!(A  h)Au„(zi)|.  Since 


is  nonsingular  if  h  j*  0,  it  follows  that  N(h) £  0. 
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i" 


l 

Suppose  that  M(zj,h)  =  0.  Taking  the  lialt  as  h  ■*  0,  |Ah  ^(z^)!  becoaes 

.  <*>  , 

Iuq  (aj)|,  the  Hronskian  of  u^,...,  Ug  evaluated  at  z  ■  z^. 

The  Hronskian  of  a  set  of  linearly  independent  analytic  functions  cannot 
vanish  Identically.  Thus,  except  for  at  dost  a  finite  set  of  values  of  z2, 

!“!<*!>  “•  “!<*<!>  I 


0 


ud(*i)  ...  ud(z<i) 


for  at  Most  a  finite  set  of  values  of  h. 

Therefore,  in  practice,  for  almost  all  pairs  (z2,h),  lu^z*)!  /  0,and  B 

2  2 

is  nonsingular.  When  B  is  nonsingular  we  have,  where  ^  -  — « - j -  , 

dzZ  c  (z) 


1  1 

...  1  ^ 

j  f<*l> 

f(*2> 

...  f(zd)  \ 

-It2  -It2 

^1  k2 
a 

e 

**•  (*k]) 

e 

e 

e 

a 

lj»f(z2) 

a 

a 

...  V»f(*d)  1 

a  a  I 

•  a  | 

• 

(-k2)d  (-kj)d 

e 

•••  <”kd)d  j 

4 

a 

1  ♦df(«1) 

a 

*df(«2) 

a  a  J 

...  ij;df(z<j)  j 

(Notice  that  both  of  the  matrices  which  have  been  written  out  have  dli 
slon  d  +  1  by  d  while  B~1  has  dimensions  d  by  d.)  This  implies  that  each 
d  +  1  by  1  column  vector 


I  <\ 
# 


can  be  written  as  a  linear  combination 


n- 


of  the  d  different  d  +  1  by  1  column 


vectors 


f(*1> 

ij)f(Zj) 


<J>df  (*  j) 


Therefore 


1  f(*l)  • 

y  • 

e  e 

e  e 

yd  <J>d  f(*!>  • 

o 

has  exactly  the  roots  y  -  -  kj. 


.  f(*d> 

.  ijjf(rd) 

e 

e 

.  *df(*d> 
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Appendix  II 

We  shall  derive  an  alternative  to  the  equation  obtained  In  Appendix  I, 
under  the  conditions  that  each 

Un(i)  a  eln  ($n  +  9nz) 

d2  -2 

for  constants  <j>n  and  0n.  As  an  approximation  to q  -  we  use  h  z  where 

-2  -2  dZ 
h  Zkg(*)  -  h  2(g(z  +  h)  -  2g(z)  +  g(z-h)). 


Then  h'2^  e 


l6n*  /  ««*>  -10nh\ 

k  e  -  h  ^e  -2  ♦  e  )e 

-  )'^v| , 


sin  “  e  h 
2  n 


(As  a  check.  If  h  ‘♦0  we  obtain  -  02  e  n  which  agrees  with  — «-  e1  nZ*) 

dzz 

Since  we  are  not  going  to  pass  to  the  limit,  we  discard  the  h*~2  factor 


and  shall  only  use  In  what  follows.  Writing 

l/i<*n+«n*>  -U*n  *  ®nf  A 
8in(0n+0n*)-  TT\«  "•  / 


8in(0n+0n*)-  TT\«  "•  / 

we  see  that 

♦l,  8ln(<J>n  +  0nz)  -  -<2  sin(<frn|h)2  sin(<|>n  +  0nz). 

.2  2 

The  exact  form  of  the  operation  <!'■'  — ~  -  -•=■ -  used  In  Appendix  I  was 

dz^  c  (z) 

irrelevant  to  the  argument  there  leading  to  the  polynomial  equation  for  the 

-k2.  What  was  Important  was  that  the  operator  used  had  the  un(z)  as 

eigenfunctions  (and  that  the  eigenvalues  were  distinct).  Thus  we  conclude 

1  9 

that  if  h  is  chosen  so  that  the  (ain©njh))£  are  distinct,  then 
|l  f(z»)  ...  f(»A)  I 


f(*l> 

•••  f(*d) 

♦h*(*l) 

a 

...  lj>nf(zd) 

a 

e 

• 

•••  ♦hf<*d> 

1  9  1 

has  as  its  roots  exactly  the  -(2  aln|0_  =  h)  . 


Appendix  III 

As  before  It  Is  assuaed  that  f(z)  has  been  aeasured  at  equally  spaced 
points  zi,...,zn  with  zj+i  -  zj  ■  h.  Next  below  define  f(s)  to  vanish 
Identically  If  z  <  zj. 

Let  r  >  1.  For  1c  -  1,2,...,  and  1  <  1  <  N,  set 


■>»»  •  ,L  (”!'■) 


f (zl-k-j) * 


Notice  that 


-U 


VW-  r  Fk<z  >  - 


+  f  (z 


i+l-k) 


-  (.f0  (kj+lj)  f(*  -(k-l)-(j+l)>r"j"^+  f(*A+l-k>» 


which  equals  f(z^)  if  k  -  1,  but  which  otherwise  Is  equal  to 


(k  L)f  <*l-(k-l)-J)r”J  *  Fk-l(*t># 


Setting  F0(za)  -  f(z£)  we  have  that  if  k  )►  1,  Fk(zJl+1)-r"1Fk(z]l)  -  F^C^) 
The  Idea  Is  to  choose  r  sufficiently  large  that  r“J  is  quite 

j 

small  whenever  A  -  k  -  J  £  0.  Then  each  Ffc(z  )  will  "almost”  be  the  same 
as  the  value  which  would  have  been  obtained  from  using  the  (unknown)  repre¬ 
sentation  of  f(z)  as  £  A,,  sin  ($n  +8n*)  to  define  f(z)  for  z  •  z0,z_^,*** 

n 

Notice  that  00 

£ 

j“0 


io  m  r’3e  i(*» + vww' 

(j"o(k+^1)rJe"i<en<J+,')h,)ei<VV) 


.  ie  h  .  .  i^n+0nZ) 

rk(re  n  -1)  k  e  n  n 
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Setting  each  rCre*^*11-!)'1  -  Rn  exp(i*n)>  every 

W  =  &  An  RS  flln  <*n  +  k  *n  +  6n  *A>‘ 

At  this  stage  not  knowing  the  8n's  and  $a's  [which  are  what  we  wish  to 

determine]  all  we  can  really  use  of  the  above  formula  is  that  each  FfcCz^)  la 
a  linear  combination  of  2d  exponentials. 

By  an  argument  quite  similar  to  that  in  Appendix  I,  but  with  the  2d 
exponentials  here  replacing  the  d  sinusoidal  functions  there,  we  see  that 
the  roots  of  the  polynomial  equation 

1  Fo(2$,)  •**  F2d-l(z8) 

y  Fi  (z£>  •••  F2d(*i,) 


-  0 


are  y 


y2d  F2d<zJl> 


+ie  h  . 

r(re~  n  -l)"1  -  R^1. 


F4d-1^*  2? 
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Appendix  IV 

There  is  a  mathematical  theory  of  best  rational  (function)  approximation 
to  power  series.  Assuming  a  property  of  g(x)  (the  power  series  to  be  approxi¬ 
mated)  which  is  called  normality  (because  it  "normally"  holds)  the  following 
are  true: 

For  each  pair  of  non-negative  integers  m  and  n  there  exist  unique  monic 

polynomials  pmn(x)  and  qmn(x)  having  degrees  exactly  m  and  n,  respectively, 

such  that  Pnm(*)(qmn(K))’1  agrees  with  g(x)  *■  E  a.  x  in  its  power 

J»0  J 

series  expansion  about  x  »  0  through  the  term  a  .  x®*11  but  no  further. 


M 

If  we  consider  a  polynomial  gM(x)  ■  £  x^  its  coefficients  for 

J  >  M  are  not  typical  of  the  class  of  all  power  series  in  x.  We  do  not 
expect  g^x)  to  be  normal;  however,  one  expects  a  random  g^(x)  to  display 
the  attributes  of  a  normal  series,  as  long  as  the  coefficients  of  powers 
of  x  higher  than  M  are  not  involved. 

The  pmn  q^'1  are  known  as  the  Pade'  approximation  (of  g(x)  or  gN(x))  and 
they  are  denoted  by  rmn,  l.e.  rB  •  p^j  q^*1.  In  1966  Wynn  discovered  the 
following  new  identify  relating  the  r^  for  a  normal  g(x).  Wynn  showed  that 

^rm+l,n  “  rm,n^  +  ^rm-l,n  ”rm,n^  *  ^rm,n+l  “rm,n^  +  ^rm,n-l  ~ 


v-1 


m.n-’ 


We  quote  from  [1]: 

"Wynn's  identity  permits  the  recursive 
table  from  the  partial  sums  r^x)  -  Ca(x) 
scheme  in  Table  7. 


construction  of  a  normal  Fade' 
■  . 

■  £  ojx-J  according  to  the 


J-0 
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Table  7 


(Notice  that  to  Initialise  properly  each  rBj„i  *  and  each  r_i>n  -  0. 
The  calculation  proceeds  in  the  order  Indicated  by  the  arrows.) 

Actually  calculating  these  ra>n  as  functions  night  be  very  difficult 
unless  one  possessed  syabol  manipulation  software.  On  the  other  hand  upon 
replacing  x  by  any  specific  value  x0  in  the  r->0  »  Ca  the  calculation  yields 

rm,n  <*o>- 


Appendix  V 

Ve  assume  that  the  u^z)  which  appear  in  the  modal  decomposition  of  f(z) 
are  known  and  that  each  modal  satisfies  a  linear  differential  equation  of  the 
form 


/sL  .  _sL_  +  k2  \  „ 

\dz2  c2(z>  "  )  " 


(z)  -  o. 


Since  the  modals  are  known  the  set  of  constants  kn  is  known;  what  is  unknown 
Is  the  set  of  coefficients  An  of  the  un(z)  in  the  representation  of  f(z). 

Suppose  for  simplicity  that  modals  u^  ***,  u^  appear  in  the  representation 


of  f.  Then  one  notes  that 


l  (kf  -  kj)  An  u„  (z) 
n-1 


d 

I 

n-2 


-  I  (k|  -  kj)  A„  un(z). 


Applying  - f -  +  k^  to  gj/z)  yields 

\dzz  C  (z)  y 

<Z)  'j3  (“2  -  “n )  (k2-kn)  A»U»W- 


g2 

Continuing, 


g. 


«  •[%({$  ■  fc  +k")] 

■[jj  (kn  -  kd)]Ad  udW- 


f(z) 


90 


n~  * 


Generalizing, 


n  (k2  -  k2\  A  u.(z)  -  II  | 

n/J  3  J  j*>\ 


.2  2 

d  a) 


.  2  2,  . 

dz  c  (z) 


0 


+  K  f(z>. 


for  n-  l,2,***,d. 

One  night  well  worry  about  calculating  derivatives  of  f  numerically  up 
to  the  order  2d~2,  if  d  is  any  size.  One  possible  way  out  is  the  following 
approach: 

Set  x  ”  tCz-z^Xzjj-z^)”^.  Notice  that  the  equation 

8i 


!i  ■  (i  -  i  +  k0 £<I) 


implies,  by  multiple  use  of  integration  by  parts,  that 


Tn 

1  «!<*> 


sin  x  dz 


sin  x  dz 


t  (^)H 


sin  x  dz  + 


f. 


N  f(l)  d  .  2  .  . 
f  (z)  d^  (8in  x)  dx 


/  f(z>  [4-  +  ) 

\  \c(z)  1  ) 

f<i>fc+k0  -1-* 


7. 

L 


N  2 

f(z)  — o  (sin  x)  dx 
dz 


sin  x  dx 


t(,)  (A  -  ~=r~  +k0 

•'•i  \d*  c(z)  1/ 


sin  x  dx. 


SI 


One  see  that  In,  general.  If  we  have  d  modal s  In  the  decomposition  of 


f(z)  then 


&-7i+k0]fW) 


2d-2 

sin  x  dz 


2d-2 

sin  x  dz 


■  A"/"  JJ,(kJ“2Kw 


2d-2 

sin  x  dz. 


One  hopes  that 


r 

J  un(z) 


2d-2 

sin  x  dz 


is  not  too  small  and  that  An  can  be  determined  upon  using  numerical 
integration  to  approximate 

S 

i-  i  7  7  v  -i 

1-2. 


rN 

J  f(I)rn  .  j  +k2\isl^ 

Z1  [j^n^dz  c  (z)  J/J 


dz 


from  knowledge  of  the  value  of  the  integrand  at  the  points  Zj,  Z2»*’*» 

With  the  basic  idea  of  the  above  derivation  in  mind,  it  becomes  clear 

2d-2 

that  we  did  not  have  to  choose  sin  x  above  as  a  multiplier.  Other  possible 
2d-2  2d-2+k 

choices  are  sin  x  un(z),  sin  x  for  a  positive  integer  k,  or 
sin  x  un(x),  also  for  a  positive  integer  k.  We  refer  to  these 

possibilities  as  option  one,  option  two,  and  option  three. 

Option  one  has  the  good  feature  that  now  the  coefficient  of  An  is  an 
integral  having  an  integrand  which  always  maintains  the  same  sign.  Option 
one  has  the  potential  drawback  that  if  lun(z)|  <  1  and  u^z)  does  not 
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change  sign  on  [zi,Zty]  we  have  not  improved  our  situation.  A  further 
drawback  is  that  calculating 


r.”  (^-4-  +k2V 

LJ'n  \dz2  c2(z)  3/j 


2d-2 

sin  x  u  (z) 
n 


may  be  much  more  difficult  than  calculating  the  same  expression  without  the 
un(z). 

Call  the  integrand  obtained  after  many  integrations  by  parts  G^(z). 

Option  two  forces  G^-^  =  for  j  “  0,1,...,  k-1.  For  such  integrands 

N-l 

the  very  simple  Integration  formula  Gk(zi)  can  be  expected  to  be  quite 

j-1 

good.  (One  can  give  an  argument  as  to  why  this  is  so  which  is  based  on  the 
rate  of  convergence  of  the  Fourier  series  for  Gk(z)  on  [zj,  z^j] .  The  co¬ 
efficients  of  e+2inx  should  be  smaller  in  absolute  value  than  a  constant, 
depending  on  k,  times  n-^.)  This  is  a  good  feature  of  option  two;  the  only 
obvious  bad  feature  is  that  the  integrand  gets  somewhat  more  complicated  as 
k  increases.  Option  three  simply  combines  options  one  and  two. 

Suppose  that  c(z)  is  slowly  varying  enough  to  be  well  approximated  on 
[  z^,zN]  by  a  constant  and  that  k  is  even,  say  k  *  2k^.  Then  (sin^x)^'"1+k1  = 

(1  cos  2  x  \  d-1+ki 

2 - ~2 — ~  )  »  and  appears  that  the  problem  of  evaluating 


r  n  /dj _ ojL  +  k2V 

c2(z)  VJ 


(sin^x)^  z  „  8  Z2»*”>  zN  can  be  carried 


out  inexpensively  using  the  Fast  Fourier  transform  (FFT)  on  [z]i,z^}. 
What  about  options  one  and  three  where  one  must  evaluate 


.*■"  (l2'c“(z)+k0. 


.  2d-2+2k,  ,  . 

sin  x  1  un(x). 


for  k,  >  0,  at  z  ■  z^,  Z2,’*’,  zN?  This  could  be  difficult  without  either 
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software  capable  of  symbol  manipulation  or  more  mathematical  investigation. 

One  possible  trick  is  the  following:  suppose  positive  Integers  and  J2  can 

be  found  with  only  a  little  larger  than  one  and  j2  only  a  little  smaller 

than  N  such  that  un(z)  ■  sin(@nz)  is  periodic  with  a  period  of  (approximately) 

zd  -  .  Then,  substituting  the  interval  [zdl»zd2l  f°r  [z^msg]  above 

2  1 

(Including  in  the  definition  of  x),  we  have  that  un(z)  ■  sin  2^x  for  some 
positive  integer  !■ .  In  this  case  the  FFT  (but  now  the  FFT  on  [z^  ,  z<j  ]) 


can  be  used  to  calculate 


1  2 


2d-2+2k 

x 


1  u  (x)  also, 
n 
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